Power Transfer in Non-linear Gravitational Clustering and Asymptotic Universality 
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We study the non-linear gravitational clustering of coUisionless particles in an expanding back- 
ground using an integro-difFerential equation for the gravitational potential. In particular, we address 
the question of how the non-linear mode-mode coupling transfers power from one scale to another 
in the Fourier space if the initial power spectrum is sharply peaked at a given scale. We show that 
the dynamical equation allows self similar evolution for the gravitational potential </f>k(i) in Fourier 
space of the form <^k(i) ~ Fi't)D(k) where the function F{t) satisfies a second order non-linear dif- 
ferential equation. We provide a complete analysis of the relevant solutions of this equation thereby 
determining the asymptotic time evolution of the gravitational potential and density contrast. The 
analysis shows that both F{t) and £)(k) have well-defined asymptotic forms indicating that the 
power transfer leads to a universal power spectrum at late times. The analytic results are compared 
with numerical simulations showing good agreement. 



If power is injected at some scale L into an ordinary 
viscous fluid, it cascades down to smaller scales because 
of the non-linear coupling between different modes. The 
resulting power spectrum, for a wide range of scales, is 
well approximated by the Kolmogorov spectrum which 
plays a key useful role in the study of fluid turbulence. 
It is possible to obtain the form of this spectrum from 
fairly simple and general considerations though the ac- 
tual equations of fluid turbulence are intractably compli- 
cated. 

In this paper, we address the corresponding question 
for non- linear gravitational clustering of coUisionless par- 
ticles in an expanding background: If power is injected 
at a given length scale very early on, how does the dy- 
namical evolution transfer power to other scales at late 
times? In particular, does the non-linear evolution lead 
to an analogue of Kolmogorov spectrum with some level 
of universality, in the case of gravitational interactions? 

We will show that the answer is essentially in the af- 
firmative. If power is injected at a given scale L = 27r/fco 
then the gravitational clustering transfers the power to 
both larger and smaller spatial scales. At large spatial 
scales the power spectrum goes as P{k) oc k'^ as soon as 
non-linear coupling becomes important. (This result is 
known in literature; see for example 0,0] ■) More inter- 
estingly, the cascading of power to smaller scales leads 
to a universal pattern at late times just as in the case 
of fluid turbulence. By studying the relevant equations 
we will show that the gravitational potential has the form 
4''k{t) = F{t)D(k) where F{t) satisfies a non-linear differ- 
ential equation and -D(k) satisfies an integral equation. 
We analyze the relevant equations analytically as well 
as verify the conclusions by numerical simulations. The 
study confirms that non-linear gravitational clustering 
does lead to a universal power spectrum at late times if 
the power is injected at a given scale initially. 

While the gravitational evolution of coUisionless par- 



ticles in an expanding background is of considerable im- 
portance in cosmology, it is probably fair to say that our 
analytic understanding of this problem is rather patchy. 
(For a general review of statistical mechanics of gravitat- 
ing systems, see previous work on analytic approx- 
imations include Zeldovich-like approximation per- 
turbative techniques 's^ , and non- linear scaling relations 
^] among many other approaches.) In cosmology there 
is very little motivation to study the transfer of power 
by itself and most of the numerical simulations in the 
past concentrated on evolving broadband initial power 
spectrum. 

We shall now provide the mathematical details of our 
formalism. If Xi(t) is the trajectory of the i—th particle, 
then equations for gravitational clustering in an expand- 
ing universe in the Newtonian limit can be summarized 
as 
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where Pb{t) is the smooth background density of matter. 
Usually one is interested in the evolution of the density 
contrast <5(t, x) = [p(i,x) — Pb{t)]/ Pb{t) rather than in 
the trajectories. Since the density contrast in the Fourier 
space Sk{t) can be related to the trajectories of the parti- 
cles, one can obtain an equation for 6k{t) from the above 
equation. This equation, in turn, will involve the veloc- 
ities of the particles and hence will not be closed. It is 
however possible to provide a closure condition for this 
equation using the following two facts, (a) At any given 
time t, particles which are already part of a bound, viri- 
alized cluster do not contribute significantly to the non- 
linear terms, (b) The velocities of the remaining particles 
can be very well approximated by the Zeldovich ansatz 
which takes the velocities to be proportional to the gra- 
dient of the gravitational potential. Given this ansatz, 
it is possible to write down a closed intcgro-differential 
equation for the evolution of gravitational potential 4'k(t) 
in the Fourier space. (The details of this derivation can 
be found in the companion paper and will not be re- 
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FIG. 1: The solution to Eq. is plotted in both the g — a plane and the u — g plane. The function g{a) asymptotically 
approaches unity with oscillations which are represented by the spiral in the right panel. The different curves in the left panel 
correspond to the rescaling freedom in the initial conditions. The entire family of solutions is represented by a single curve in 
the u — g plane. One fiducial curve which was used to model the simulation is shown in the red. 



peated here.) The final equation is 
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This equation governs the dynamical evolution of the 
system but, of course, is intractably complicated. In ad- 
dition to the obvious difficulty of solving an integro dif- 
ferential equation, we also need to tackle the issue of 
incorporating appropriate initial conditions, which will 
influence the evolution at the early stages. (Some re- 
sults based on perturbation solution to this equation are 
discussed 0-) Our aim is to look for late time scale 
free evolution of the system exploiting the fact that the 
above equation allows self similar solutions of the form 
(j)]^{t) = F(i)L)(k). Substituting this ansatz into Eq. Q 
we obtain two separate equations for F{t) and -D(k). It 
is also convenient at this stage to use the expansion fac- 
tor ait) = {t/to)'^/'^ of the matter dominated universe as 
the independent variable rather than the cosmic time t. 
Then simple algebra shows that the governing equations 
are 
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Equation Q governs the time evolution while Eq. Q 
governs the shape of the power spectrum. (The sepa- 
ration ansatz, of course, has the scaling freedom F — > 
fj,F,D — > (l//i)Z? which will change the right hand side 
of Eq. © to -/iF^ and the left hand side of Eq. Q to 
fiHQDk. But, as to be expected, our results will be inde- 
pendent of fi; so we have set it to unity). Our interest lies 
in analyzing the solutions of Eq. ^ subject to the initial 
conditions F{ai) — constant, (dF/da)i = at some small 
enough a — ai. 

Inspection shows that Eq. (j2Jl has the exact solution 
F{a) = (3/2)a~^. This, of course, is a special solution 
and will not satisfy the relevant initial conditions. How- 
ever, Eq. lI2Jl fortunately belongs to a class of non-linear 
equations which can be mapped to a homologous system. 
In such cases, the special power law solutions will arise 
as the asymptotic limit. (The example well known to as- 
tronomers is that of isothermal sphere J]. Our analysis 
below has a close parallel.) To find the general behaviour 
of the solutions to Eq. Q , we will make the substitution 
F{a) = {3/2)a~^g{a) and change the independent vari- 
able from a to q — log a. Then Eq. ^ reduces to the 
form 
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This represents a particle moving in a potential V{g) — 
{l/2)g^ — (3/4)g^ under friction. For our initial condi- 
tions the motion will lead the "particle" to asymptoti- 
cally come to rest at the stable minimum at g = 1 with 
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FIG. 2: Left panel: The results of the numerical simulation with an initial power spectrum which is a Gaussian peaked at 
L = 24. The y-axis gives Afc where A| = k^P{k)/2n'^ is the power per logarithmic band. The evolution generates a well known 
tail at large scales (see for example, and leads to cascading of power to small scales. Right panel: The simulation data 
is re-expressed by factoring out the time evolution function g{a) obtained by integrating Eq. 0. The fact that the curves fall 
nearly on top of each other shows that the late time evolution is scale free and described by the ansatz discussed in the text. 
The rescaled spectrum is very well described by P{k) oc fc^'''*/(l + {k/kof"^) which is shown by the completely overlapping, 
broken blue curve. 



damped oscillations. In other words, F{a) (3/2)a~^ 
for large a showing this is indeed the asymptotic solu- 
tion. From the Poisson equation in Eq.QJ, it follows 
that fc^0j^ oc so that 5-^{a) oc g{a)k'^ D{k) giv- 

ing a direct physical meaning to the function g{a). The 
asymptotic limit corresponds to a rather trivial case of 
5-^ becoming independent of time. What will be more in- 
teresting — and accessible in simulations — will be the 
approach to this asymptotic solution. To obtain this, we 
introduce the variable v = 2{dg/dq) so that our system 
reduces to the "phase space" equations 

1 V 

^^-^v-Sgig-'^y, .9=2 (6) 

where the dot denotes differentiation with respect to q. 
Dividing the first equation by the second and changing 
variables to u ~ {v + g), we get the first order form of 
the autonomous system to be 

^ ^ 6g(.9 - 1) 
dg g 

The critical points of the system are at (0, 0) and (1, 1). 
Standard analysis shows that: (i) the first one is an un- 
stable critical point and the second one is the stable criti- 
cal point; (ii) for our initial conditions the solution spirals 
around the stable critical point. 

Figures 1(a) and 1(b) describe the solution in the g — a 
and u — g planes. The g(a) curves clearly approach the 



asymptotic value of g « 1 with superposed oscillations. 
The different curves in Fig 1(a) are for different initial 
values which arise from the scaling freedom mentioned 
earlier. (The thick red line correspond to the initial con- 
ditions used in the simulations described below.). Fig 
1(b) shows this family of solutions in the u — g plane. As 
usual, in going from the second order equation for g to 
the first order equation in u — g plane, we map a fam- 
ily of solutions to a single curve The stable critical 
point acts as an attractor. The solution g{a) describes 
the time evolution and solves the problem of determining 
asymptotic time evolution. 

To test the correctness of these conclusions, we per- 
formed a high resolution simulation using the TreePM 
method 0,13 and its parallel version ^3 with 128"^ par- 
ticles on a 128'^ grid. We used a cubic sphne softened 
force with softening length e = 0.4 to ensure collisionless 
evolution in the simulations. We computed all relevant 
statistics at scales L larger than 2e to avoid errors due 
to the force softening. Details about the code param- 
eters can be found in Q. The initial power spectrum 
P{k) was chosen to be a Gaussian peaked at the scale 
of hp = 271 /Lp with Lp = 24 grid lengths and with a 
standard deviation Afc = 27r/Lf,oa;, where Li^ox = 128 
grid lengths is the size of one side of the simulation vol- 
ume. The amplitude of the peak was taken such that 
Ah„ (fc = A:p,a = 0.25) = l. 
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The late time evolution of the power spectrum (in 
terms of A^. = k^P{k)/2Tr^ where P = jJ^p is the power 
spectrum of density fluctuations) obtained from the simu- 
lation is shown in Fig 2(a). In Fig 2(b), we have rescaled 
the Afe, using the appropriate solution g{a). The fact 
that the curves fall on top of each other shows that the 
late time evolution indeed scales as g{a) within numerical 
accuracy. A reasonably accurate fit for g{a) used in this 
figure at late times is given by 



This index has a simple interpretation [ill IT^ . The 
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The key point to note is that the asymptotic time evo- 
lution is essentially 5{a) oc a except for a logarithmic 
correction, even in highly non-linear scales. (This was 
first noticed from somewhat lower resolution simulations 
in jl]|.). Since the evolution at linear scales is always 
(5 oc a, this allows for a form invariant evolution of power 
spectrum at all scales. Gravitational clustering evolves 
towards this asymptotic state. 

This behaviour is fairly generic and we have performed 
a series of simulations with different initial conditions to 
test this claim. Fig. 3 shows scaled by g{a) for a 
simulation where we have initial power concentrated in 
two narrow windows in fc-space. In addition to power 
around Lp — 24 grid lengths as in the previous case, we 
added power at fci — 2tt/Li {Li — 8 grid lengths) using 
a Gaussian with the same width as before. Amplitude 
at Li was chosen to be five times higher than that at 
Lp, which means that Aii„{ki,a — 0.05) = 1. The fact 
that the curves fall on top of each other in this case too 
shows that the late time evolution indeed scales as g{a), 
independent of initial conditions. 

The form of the power spectrum in Fig. 2 is well ap- 
proximated by 
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This fit is shown by the broken blue line in the figure 
which completely overlaps with the data and is barely 
visible. (Note that this fit is applicable only at L < Lp 
since the k'^ tail will dominate scales to the right of the 
initial peak; see the discussion in (111]). At non-linear 
scales P{k) oc k~^ making Afe fiat, as seen in Fig. 2. (This 
is not a numerical artifact and we have sufficient dynamic 
range in the simulation to ascertain this.) At quasilinear 
scales P{k) oc The power spectrum in Fig. 3 is 

also well-approximated by Eq. excepting for the fact 
that the scale 27r/fco ~ 4.1 in this case. A lower value for 
the scale is to be expected because of the additionally 
injected power at the lower scale Li. 

The effective index of the power spectrum varies be- 
tween —3 and —0.4 in this range of scales for which the 
fitting function has been provided. To the lowest order of 
accuracy, the power spectrum at this range of scales is ap- 
proximated by the mean index n « —1 with P{k) oc k~^. 




100 



L=2 7r/k 



FIG. 3: The results of a numerical simulation with an initial 
power spectrum which has two Gaussians peaked at L — 24 
and L = 8. The figure has been plotted by factoring out 
the time evolution function g{a) (Eq. 0) from A^. The fact 
that the curves fall nearly on top of each other in this case 
as well shows that the late time evolution is scale-free and 
is described by the ansatz discussed in the text, independent 
of initial conditions. The rescaled spectrum is once again 
described by P{k) oc fc~'''^/(l + (k/ko)^'^) which is shown by 
the completely overlapping, broken blue curve. 



ensemble-averaged gravitational potential energy of fiuc- 
tuations per unit volume is given by 
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This energy reaches equipartition (i.e., contributes same 
amount per logarithmic band of scales) when P{k) oc 
k^^. The same result holds for kinetic energy if the mo- 
tion is dominated by scale invariant radial fiows. Our 
result suggests that gravitational power transfer evolves 
towards this equipartition. In principle, this k depen- 
dence of the power spectrum is determined by Eq. I^J ; but 
analytical understanding of this equation is more difficult 
since the solution depends on the phases of which do 
not contribute to the power spectrum. This issue is under 
investigation. 
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